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Abstract 

In this paper we present a classical Monte Carlo simulation of the orthorhom- 
bic phase of crystalline polyethylene, using an explicit atom force field with 
unconstrained bond lengths and angles and periodic boundary conditions. We 
used a recently developed algorithm which apart from standard Metropolis lo- 
cal moves employs also global moves consisting of displacements of the center 
of mass of the whole chains in all three spatial directions as well as rotations 
of the chains around an axis parallel to the crystallographic c-direction. Our 
simulations are performed in the NpT ensemble, at zero pressure, and ex- 
tend over the whole range of temperatures in which the orthorhombic phase 
is experimentally known to be stable (10 - 450 K). In order to investigate 
the finite-size effects in this extremely anisotropic crystal, we used different 
system sizes and different chain lengths, ranging from C12 to Cge chains, the 
total number of atoms in the super-cell being between 432 and 3456. We 
show here the results for structural parameters, such as the orthorhombic 
cell parameters a, b, c, and the setting angle of the chains, as well as internal 
parameters of the chains, such as the bond lengths and angles. Among ther- 
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modynamic quantities, we present results for thermal expansion coefficients, 
elastic constants and specific heat. We discuss the temperature dependence 
of the measured quantities as well as the related finite-size effects. In case of 
lattice parameters and thermal expansion coefficients, we compare our results 
to those obtained from other theoretical approaches as well as to some avail- 
able experimental data. We also suggest some possible ways of extending this 
study. 
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I. INTRODUCTION 



Macromolecular crystals represent at the same time a particular case of a solid state 
system as well as that of a polymer system. The simplest and paradigmatic case in this 
field is crystalline polyethylene. As is well known, it is very difficult to obtain reliable 
experimental data for these systems, mainly due to problems associated with preparing 
sufficiently large mono crystalline samples. This increases the value of theoretical insight for 
understanding and eventual prediction of their properties. However, compared to the case 
of liquid polymer phases, a quantitative description of structure and properties of a solid 
phase is a more subtle problem. While in the former case it is often possible to replace whole 
monomer groups by united atoms and constrain the bond lengths and sometimes even the 
bond angles, in case of macromolecular crystals, instead, it is preferred to employ a force 
field for all atoms, without enforcing any constraints on the degrees of freedom [Q]. 

Apart from the study of the ground state crystal structure and properties, which can be 
readily done within the framework of molecular mechanics, once a force field is available, 
the main interest lies in the finite temperature properties. At low temperatures, where 
the creation of conformational defects is unlikely, and the displacements of the atoms are 
relatively small, the structure of the PE crystal is dominantly governed by packing energetics 
of all-trans conformation chains. In this regime, phonon modes are well defined excitations of 
the system, and a possible theoretical approach is the use of quasi-harmonic or self-consistent 
quasi-harmonic approximations, which have the advantage of allowing for quantum effects 
to be taken into account easily [H-^j. 

Such methods are, however, intrinsically incapable of treating correctly the large am- 
plitude, anharmonic motion of whole chain segments, which arises as the temperature ap- 
proaches the melting point (at normal pressure about 414 K for a PE crystal). Actually, for 
crystalline PE they start to fail just at the room temperature of 300 K If the melting is 
prevented by increasing the external pressure, the orthorhombic crystal structure may per- 
sist to higher temperatures, and eventually undergo a phase transition into the hexagonal 
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"condis" phase, characterized by a large population of gauche defects 0. In order to study 
the complex behavior in this high temperature regime, computer simulation techniques, 
like molecular dynamics or Monte Carlo, appear as a particularly convenient tool. Using a 
suitable ensemble and an appropriate simulation technique one can evaluate various ther- 
modynamic quantities, like specific heat, elastic constants or thermal expansion coefficients, 
and structural quantities, such as bond lengths and angles, or defect concentration. It is 
also possible to directly study a variety of physically interesting phenomena associated with 
the structural phase transitions between different crystal modifications. 

Computer simulations are, however, still limited to relatively small systems. In case 
of macromolecular crystals, the situation is rather peculiar due to the extreme anisotropy 
of the system originating from its quasi-one-dimensional character. When the chains are 
short, the system crosses over from a PE crystal to an n-paraffin crystal, which behaves 
in a substantially different way. Due to an easy activation of chain rotation and diffusion, 
ra-paraffin crystals undergo, as the temperature is increased, a characteristic series of phase 
transitions, depending on the chain length ("rotator phases", This fact makes it nec- 

essary to study and understand the related finite-size effects, if the results of the simulation 
are to be regarded as representative for the limit of very long chains, corresponding to PE, 
and interpreted in a consistent way. The choice of boundary conditions also represents a 
non-trivial issue, as documented by the very extensive explicit atom MD simulation in 
which a transition from an initial orthorhombic structure to a parallel zig-zag chain arrange- 
ment was observed already at 111 K, contrary to experiment, probably because of the use 
of free boundary conditions in all spatial directions and a related large surface effect. 

Finite temperature simulations of realistic, explicit atom models of crystals with long 
methylene chains, of which we are aware so far, have exclusively used molecular dynamics 
as the sampling method. Ryckaert and Klein used a version of the Parrinello-Rahman 
MD technique with variable cell shape to simulate an ra-paraffin crystal with constrained 
bond lengths and periodic boundary conditions. Sumpter, Noid, Liang and Wunderlich 
have simulated large crystalhtes consisting of up to 10^ CII2 groups, using free boundary 
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conditions [Q]. Recently, Gusev, Zehnder and Suter ^ have simulated PE crystal using 
periodic boundary conditions, at zero pressure, using also the Parrinello-Rahman technique. 

On the other hand, not much MC work seems to exist in this field, and to our knowledge, 
the only work is that of Yamamoto which was aimed specifically at the high tempera- 
ture "rotator" phases of n-paraffins. It concentrated exclusively on the degrees of freedom 
associated with the packing of whole chains, assumed to be rigid, neglecting completely the 
internal, intramolecular degrees of freedom. 

The aim of this paper is twofold. On the one hand, we would like to explore the applica- 
bility of the MC method to a classical simulation of a realistic model of a PE crystal, using 
an explicit atom force field without any constraints, and periodic boundary conditions in 
all spatial directions. In order to have a direct access to quantities like thermal expansion 
coefficients, or elastic constants, we choose to work at constant pressure. If MC turns out 
to be a well applicable method for this system on the classical level, it would be a promising 
sign for an eventual Path Integral Monte Carlo (PIMC) study allowing to take into account 
also the quantum effects, known to play a crucial role at low temperatures 0-^]. It is known 
that for path integral techniques the use of MC is generally preferred to MD, because of the 
problem of non-ergodicity of the pseudo-classical system representing the quantum one in a 
path integral scheme. As a related problem, we would also like to understand the finite-size 
effects and determine for each particular temperature the minimal size of the system that 
has to be simulated in order to be reasonably representative of the classical PE crystal. 
Such information would also be very useful for an eventual PIMC study, where an addi- 
tional finite-size scaling has to be performed, because of the finite Trotter number. On the 
other hand, we would like to study the physics of the orthorhombic phase of PE crystal in 
the classical limit over the whole temperature range of its experimentally known existence 
at normal pressure. We stress that our goal here is not to tune the force field in order 
to improve the agreement between the simulation and experiment. Our emphasis rather 
lies on providing results calculated for a given force field with an essentially exact classical 
technique. These might in future serve as a basis for estimation of the true importance of 
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quantum effects, when these can be taken into account by a PIMC technique, as well as 
for an assessment of the range of validity of different classical and quantum approximation 
schemes. Some preliminary simulation results of this study have been presented briefly in 



recent conference proceedings ||1CI|| . 



The paper is organized as follows. In Sect. II, we describe the force field as well as the 
constant pressure simulation method used. The MC algorithm itself will be addressed only 
briefly, since it has already been discussed in detail in its constant volume version in Refs. 
||TT| , [T0|1 . In Sect. Ill, we present the results for the structural and thermodynamic properties 
of the orthorhombic PE crystal obtained from zero pressure simulation in the temperature 
range 10 - 450 K, for different chain lengths and system sizes. We discuss the temperature 
dependence of the measured quantities as well as the related finite-size effects, and compare 
our results to those obtained from other theoretical approaches, as well as to some available 
experimental data. In the final Sect. IV we then draw conclusions and suggest some possible 
further directions. 

II. CONSTANT PRESSURE SIMULATION METHOD 

In this section we describe some details of the simulation method, as well as of the force 
field used. Before doing so, we recall here the structure of the orthorhombic PE crystal 



The unit cell contains two chains, each consisting of 2 CII2 groups, giving a total of 
12 atoms per unit cell. The all-trans chains extend along the crystallographic c-direction 
(z-axis) and are packed in a "herringbone" arrangement, characterized by the setting angle 
ip (angle between the xz plane and the plane containing the carbon backbone of a chain) 
alternating from one raw of chains to another between the values The packing is 

completely determined by specifying the 3 lattice parameters a, b, c as well as the value of 
To specify the internal structure of the chains, three additional parameters are needed, 
which may be taken to be the bond lengths rcn and rcc and the angle Orch- 

We have simulated a super-cell containing i x j x k unit cells of the crystal, k being 
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integers. Periodic boundary conditions were applied in all three spatial directions, in order 
to avoid surface effects. Our PE chains with backbones consisting of 2k carbon atoms are 
thus periodically continued beyond the simulation box and do not have any chain ends. 

As the study was aimed specifically at the orthorhombic phase of the PE crystal and 
we did not expect phase transitions into different crystal structures, we did not consider 
general fluctuations of the super-cell shape, otherwise common in the Parrinello - Rahman 
MD method. We have constrained the crystal structure angles 7 to the right angle 
value, not allowing for shear fluctuations. The volume moves employed thus consisted only 
of an anisotropic rescaling of the linear dimensions of the system by three scaling factors 
Si,S2,S3, which relate the instantaneous size of the super-cell to that of the reference one. 
The reference super-cell always corresponded to lattice parameters a = 7.25A, b = S.OOA, c = 
2.53A. The acceptance criterion for the volume moves was based on the Boltzmann factor 
(siS2S3)^e~^'^, where H = U + PV0S1S2S3, U is the potential energy of the system, p is 
the external pressure, Vq the volume of the reference super-cell, and the total number 
of atoms in the system. Throughout all simulations described in this paper, the external 
pressure was set to zero. 

We come now to the description of the potential. We have used the force field developed 
for the PE crystal by Sorensen, Liau, Kesner and Boyd fl^ , with several slight modifications 
of the bonded interaction. This consists of diagonal terms corresponding to bond stretching, 
angle bending, and torsions, as well as of off-diagonal bond-bond, bond-angle and various 
angle-angle terms. For convenience, we have changed the form of the expansion in bond 
angles, replacing the expressions {9 — Oq) in all terms by (cos 6* — cos6'o)/(— sin6'o). The 
original form of the vicinal bend-bend interaction, having two different force constants, fcr 
and kc for torsional angles close to trans and gauche minima, respectively, is useful for 
ground state studies, but not for a finite-temperature simulation, where the torsional angles 
may continuously fluctuate from one minimum to another. The form of this interaction has 
therefore been changed into the one used in Ref. , k cos (y9(cos 9i — cos 6*0) (cos 6*2 — cos 6*0) , 



where 61, 62 are bond angles, (p is a torsional angle, and k is a new force constant. The value 



of the latter constant was taken to be —kT/ sin^ for C-C-C-C torsions and 2^^/ sin^ 
for C-C-C-H torsions, in order to reproduce the curvature of the potential in the vicinity of 
the ground state equilibrium value of each torsional angle in the PE crystal. For H-C-C-H 
torsions, which are in the ground state in both trans and gauche minima, we took for k an 



approximate value oik = (^^ kckx / cos | cos vr ^ / sin^ 6*0, which guarantees that both original 
force constants k^ and kc are approximated in the vicinity of the respective minimum with 
the same error of about 8 %. At this point we comment on the torsional potential used. As 



the force field has been originally designed for ground-state studies, the explicit torsional 
potential for all torsions contains only the term cos 3^9, which yields zero energy difference 
between the trans and gauche minima in C-C-C-C torsions. The actual difference thus 
comes just from the non-bonded interaction superposed over 1-4 atoms, and is, according 



to Ref. [jT5[, about 340 K, which is a somewhat higher value than the generally accepted 
one of about 250 K. However, the use of periodic boundary conditions anyway inhibits the 
creation of conformational defects strongly, and therefore this difference is not likely to play 
an important role, at least in the temperature range studied. 

Concerning the nonbonded interaction, we used a spherical cutoff of 6 A for all pair 
interactions, which corresponded to interaction of a given atom with about 110 neighbors. 
The list of neighbors has been determined at the beginning of the simulation with respect to 
the reference structure, and kept fixed throughout the evolution of the system (topological 
interaction). The use of the topological interaction would clearly preclude the longitudinal 
diffusion of the chains by creating an artificial energy barrier, which might be unrealistic 
for a study of short alkanes. However, since our main interest lies in the limit of very long 
chains, this approximation, common in the study of crystals, is acceptable here, and brings an 
advantage of considerably speeding up the execution of the program. The reference structure 
was obtained by placing the ideal crystal structure, described at the beginning of this section, 
inside the reference box, taking for the setting angle and the internal chain parameters 
the values of = 43.0°, r^c = l.53QA,rcH = 1.09A,^hch = 107.4°, respectively. These 
values turned out to be to a good approximation close to their true average values throughout 
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the whole temperature range of the simulation, thus proving the consistency. Because of the 
relatively low cutoff used, we have added the long-range corrections to the non-bonded energy 
and diagonal components of the pressure tensor. These corrections have been calculated for 
the reference structure in the static lattice approximation, and tabulated for a suitable mesh 
of scaling factors Si, S2, S3. During the simulation, the values of the corrections corresponding 
to the instantaneous values of the scaling factors were calculated from the table by means of 
three-dimensional linear interpolation. We note here that our way of treating the non-bonded 
interaction is different from that used in |]l3l, where the interaction of a given atom with 



two neighboring shells of chains has been taken into account and no long-range corrections 
have been applied. 

The MC sampling algorithm for PE crystal was described in considerable detail in our 
previous papers |lTT| , p!0| , and here we present it just briefly. In addition to the volume 
moves, we used local moves to move the atoms and global moves to move the chains. In the 
local moves, the atoms of the crystal lattice were visited in sequential order, and different 
maximum displacements have been used for the atoms of carbon and hydrogen, reflecting the 
fact that a carbon atom has four covalent bonds while a hydrogen atom has just one. The 
typical value of the acceptance ratio for the local moves was kept close to 0.3, corresponding 
at temperature T = 100 K to isotropic maximum displacements of 0.03 A and 0.06 A for 
carbons and hydrogens, respectively. In the global moves, displacements of the center of 
mass of a whole chain along all three axes accompanied by rotation of the chain along a line 
parallel to the crystallographic c-direction and passing through the center of mass of the 
chain were attempted. Choosing the fraction of global moves to be 30 % in this study, the 
global and local moves were alternated at random, in order to satisfy the detailed balance 
condition. Once it was decided to perform a global move, all the chains of the super-cell 
were attempted to move in sequential order. For C12 chains at T = 100 K, the maximum 
(anisotropic) displacements of the chains were Axmax = ^Vmax = OAIA, Azmax = O.OSA, 
the maximum rotation angle being Aip^ax = 11°. This choice of the parameters resulted in 
an acceptance ratio of about 0.18 for the global moves. We note here that the energy change 
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associated with a rigid displacement or a rotation of a whole chain scales linearly with the 
length of the chain, and therefore for longer chains we reduced appropriately the parameters 
of the global moves in order to preserve the same acceptance ratio. No attempt to optimize 
the maximal displacements or fraction of different kinds of moves has been done in this 
study Concerning the volume moves, we attempted a change of all three scaling factors 
si,S2, S3 after each sweep over the lattice (MCS) performing local or global moves. For the 
super-cell 2x3x6 unit cells at T = 100 K, the maximum changes of the scahng factors used 
were 5si = 5s2 — 0.0045, 5s3 — 0.0009, where the different values reflected the anisotropy 
of the diagonal elastic constants 011,022,033. This choice resulted in an acceptance ratio of 
0.21 for the volume moves. 

We have simulated several different super-cell sizes. The smallest one consisted of 2 x 3 x 6 
unit cells, contained 12 C12 chains with a total of 432 atoms, and was used for simulation of 
the system at seven different temperatures ranging from 10 to 300 K. To study chain-length- 
dependent finite-size effects, at T = 100 K and T — 300 K the simulation was performed 
also on super-cells consisting of 12 C24 and C48 chains, containing respectively 2 x 3 x 12 and 
2 X 3 X 24 unit cells, while at T = 200 K only the latter system size was studied in addition 
to the smallest one. At T — 300 K, a super-cell size 4 x 6 x 12 unit cells, twice as large as 
the smallest one in each spatial direction, was also used, to study volume-related finite-size 
effects. Finally, at all four highest temperatures, T — 300, 350, 400, 450 K, a super-cell with 
the longest, Cge chains was used, consisting of 2 x 3 x 48 unit cells and containing 3456 
atoms. 

As the initial configuration for the lowest temperature simulation for a given super-cell 
size we used the corresponding reference structure. In course of the simulation, we made 
use of the final configuration of the run at a lower temperature when possible, and always 
equilibrated the system for at least 2 x 10^ MCS before averaging. Our statistics is based 
on the run length of 1.4 x 10^ MCS per data point for the smallest, 432 atom system, and 
a run length decreasing linearly with the number of atoms for the larger systems. 

We have calculated the specific heat at constant pressure Cp from the enthalpy fluctua- 
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tions. The pressure tensor was calculated by means of a standard virial expression. In order 
to check the consistency of the simulation algorithm, we also evaluated the kinetic energy 
from the corresponding virial expression. During the averaging, the accumulators for the to- 
tal energy, lattice parameters a, b, c and setting angle ip were updated after each MCS while 
those for virial and other quantities were updated only every 10 MCS. Histograms have been 
accumulated for the structural quantities, like the bond lengths and angles, torsional angles 
and the setting angle of the chains as well as the displacement of the center of mass of the 
chains along the 2-axis. The whole run was always subdivided into four batches and the 
batch subaverages were used to estimate the approximate error bars of the total averages. 

Concerning the elastic constants Cn, C22, C33, C12, C13, C23 (in the Voigt's notation), we have 
independently determined them in two different ways. Apart from the standard Parrinello- 



Rahman fluctuation formula 



Cik = -^{eiCk) A; = l,2,3, (1) 



we applied also the new fluctuation formula proposed in Ref . , in its approximate version 
suitable for small strain fluctuations 

Cik = -'^{Pien){enek)~^ , (2) 

n 

where Pi and Cj are the diagonal components of the pressure tensor and strain tensor, re- 
spectively. 



III. RESULTS AND DISCUSSION 

In this section we describe and discuss the results of the simulation. We start with a 
comment on the stability of the crystal structure. The initial structure with the characteristic 
"herringbone" arrangement of the chains was found to be stable at all temperatures and all 
system sizes except for the smallest system with C12 chains, where an occasional rotation 
of a whole chain was observed at T = 300 K. In the largest system with Cgg chains, no 
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change of structure was observed up to T = 450 K. Although the latter temperature is 
larger than the experimentally known melting temperature of PE crystal (414 K), the use of 
periodic boundary conditions inhibits the melting and allows the simulation of a superheated 
crystalline phase. On the other hand, our arrangement with constrained angles of the super- 
cell is compatible both with the orthorhombic and with the hexagonal phase, and would not 
prevent the system from entering the latter, which would occur when the ratio of the lattice 
parameters | reaches the value v^- Our observation thus agrees with the experimentally 
known stability of the orthorhombic structure up to the melting point. In order to appreciate 
the amount of disorder present in the system at T = 400 K, we show in Fig.l (a) a projection 
of the atoms on the xy plane for a typical configuration. The "herringbone" arrangement 
of the chains is still clearly visible, in spite of a well pronounced disorder. In Fig.l (b), a 
projection of the same configuration on the xz plane is shown. 

In Figs. 2, 3,4 we show the temperature dependence of the lattice parameters a, 6, c. Ex- 
trapolating these curves down to T = 0, we find the ground state values of a = 7.06A, b = 
4. 89 A, c = 2. 530 A. We note here that these values do not quite agree with those reported 



in Ref. [IT3, where the following values have been found a = 7. 05 A, b = 4. 94 A, c = 2. 544 A. 
We attribute these discrepancies to the different way of treating the non-bonded interaction 
as well as to our slight modifications of the bonded interaction. 

Before discussing the thermal expansion of the lattice parameters, we comment on the 
finite-size effects observed. A particularly pronounced one is observed on the lattice param- 
eter b, where the values for the smallest system with C12 chains are already at T = 100 K 
slightly larger than those for both systems with longer chains. The effect becomes stronger 
with increasing temperature. At T = 300 K, where we have data for four different chain 
lengths, the value of b is clearly seen to increase with decreasing chain length, most dramat- 
ically for the system with C12 chains. On the other hand, a considerably weaker finite-size 
effect is seen on lattice parameters a and c. While for the latter one it is perhaps not 
surprising, because of the large stiffness of the system in the chain direction, the distinct 
behavior of the b parameter with respect to the a parameter does not appear to be so 
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straightforward to interpret. We believe that its origin hes in the particular character of the 
chain packing, where the shortest hydrogen- hydrogen contact distance is just that along the 
crystallographic 6-direction (y-axis) This results in a stronger coupling of the lateral 
strain €2 along the 6-direction to the longitudinal displacements of whole chains. Since the 
fluctuations of these displacements are larger for systems with short chains, a particular 
finite-size effect arises. 

Concerning the thermal expansion itself, it is convenient to discuss separately the case of 
the lateral lattice parameters a, 6, and that of the axial one c. We start with the lateral ones, 
and show in Fig.5 also the temperature dependence of the aspect ratio |. Two regimes can 
be clearly distinguished here. For temperatures up to about 250 K, both lattice parameters 
expand in a roughly linear way with increasing temperature, and the ratio | raises only 
slightly from its ground state value of 1.44 (which differs substantially from the value of 

= 1.73, corresponding to hexagonal structure). It is characteristic for this regime that the 
thermal expansion arising due to lattice anharmonicities can be described within a phonon 
picture using a quasi-harmonic or self-consistent quasi-harmonic approximation For 
higher temperatures, the picture changes. While the lattice parameter a starts to increase 
faster, the expansion of the parameter b at the same time becomes more slow until it develops 
a maximum at T = 350 K, where it starts to decrease again. Such behavior of b has been 
already observed in the work 0. As a consequence, the aspect ratio | increases strongly. 
This suggests that the driving force of the change of the lateral lattice parameters in this 
regime is the approach of a phase transition to a hexagonal phase, in which each chain is 
surrounded by six chains at equal distance. We have actually continued our simulations 
to even higher temperatures, and from a limited amount of simulation performed in that 
region we have found an indication that the hexagonal phase is indeed reached in the range 
of temperatures 500 - 550 K (Fig.5). It is, however, clear that in order to obtain reliable 
results from the simulations in this high-temperature range, where the phase transition in 
a real PE crystal involves large amplitude displacements and rotations of the whole chains, 
as well as a considerable population of conformational defects, it would be necessary to 
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introduce several modifications into the simulation algorithm. We shall come back to this 
point again in the final section. 

In Fig. 6, the temperature dependence of the average setting angle {{ipl) of the chains is 
shown. It also fits well within the two regimes scenario, being rather fiat up to T = 300 
K, and then starting to decrease. This decrease could indicate an approach of the value of 

= 30°, compatible with the symmetry of the hexagonal phase. It is also interesting to 
note the pronounced finite-size effect, similar to that observed in the case of b. 

Before discussing the temperature dependence of the axial lattice parameter c, it is 
convenient to plot also the thermal expansion coefficients, defined as at = aC^dai/ dT,i = 
1, 2, 3, where ai, 02, 03 are the lattice parameters a, b, c. These have been obtained by taking 
the finite differences of the lattice parameters and are shown in Figs. 7, 8, 9 as a function 
of temperature. As their behavior trivially follows from that of the lattice parameters, 
discussed for i = 1,2 above, we just note here that all three coefficients converge at low 
temperatures to nonzero finite values, as can be expected in the classical limit. 

Concerning the behavior of c and a^, the characteristic feature is that 0:3 is negative in the 
whole temperature range and an order of magnitude smaller than ai. It is interesting here to 
compare our result for 03 to that found in Ref. M within the quasi-harmonic approximation 



for a different force field . A distinct feature of the latter result is that the classical value 
of 0^3 is considerably smaller in magnitude than the quantum mechanical value (and also 
the experimental one), and approaches zero as T — 0. Since in the classical limit there is 
no a priori reason for such behavior, it has to be regarded as accidental. According to the 
argumentation in Ref. 0, there are contributions of different sign to 0:3 from different phonon 
modes, negative from lattice modes (mainly backbone torsions), and positive from the harder 
ones. In the classical limit, all the modes contribute at all temperatures and happen to just 
cancel each other as T ^ 0. In order to estimate the amount of contribution of the torsions 
in our case, we have made use of the work 0], where a formula is derived for the axial 
thermal contraction in a simple one-chain model with only torsional degrees of freedom. In 
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the classical limit, the formula predicts c — cq = — |cosin^ ^{'Pcccc)^ where a = n — 9, 



ccc 



and 4>cccc is the fluctuation of the torsional angle fcccc from the trans minimum. In 
Fig. 10 we have plotted c against {(pcccc)^ found a very good linear dependence over 
the whole temperature range up to 450 K, where the effective torsional constant, deflned as 

Ctors = 777^ — \, undergoes a considerable softening with increasing temperature (Fig. 11). 

\fcccc' 

The proportionality constant determined from our plot was, however, larger in magnitude 
by about 50 % with respect to the above value of |co sin^ ^, valid for the simple model. The 
linear dependence suggests that with the force fleld used |T^, the negative as originates 



almost entirely from the C-C-C-C torsions, which points to a certain intrinsic difference in 



the anharmonic properties of the force flelds |13| and [14]. The pronounced increase of 



for T > 300 K is thus a consequence of a strong softening of the torsional potential due to 
the large amphtude of the fluctuations. 

For comparison of these results to experimental ones, we have chosen two sets of data. 



In the work [0, lattice parameters a, b, c have been measured in the temperature range 93 
- 333 K, and the data are smooth enough to allow a direct extraction of thermal expan- 
sion coefficients by means of finite differences. The other chosen set of data is to our 
knowledge the only one covering the range from helium temperatures up to T = 350 K, 
however, the scatter of the data is too large for a direct numerical differentiation. Therefore 
we decided to first fit them by a fourth-order polynomial and then evaluate the expansion 
coefficients analytically. For both sets, the experimental values of the lattice parameters are 
shown in Figs. 2, 3,4 and the corresponding thermal expansion coefficients in Figs. 7, 8, 9 . 
Concerning the absolute value of the lattice parameters, in case of a our results agree 



well with the data [|T9|, while falling slightly below the data |2^, in particular at the lowest 
temperatures. In case of b our results fall slightly above and in that of c slightly below 
both data sets in the whole temperature range. As far as the temperature dependence 
itself is concerned, this is most conveniently discussed in terms of the thermal expansion 
coefficients a^. We note first that for temperatures lower than about 150 K, quantum effects 
become crucial and cannot be neglected, as they are responsible for the vanishing of all 
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expansion coefficients in the limit T — 0. A meaningful comparison of our classical results 
to experimental data is thus possible only for larger temperatures. For T > 150 K, ai agrees 
qualitatively well with the data |T9[, although the experimental ones appear to be slightly 
smaller. In particular, the pronounced increase of ai for T > 250 K appears to be well 
reproduced, in contrast to the result obtained in Ref. within the classical quasi-harmonic 



approximation for the force field flj]. The agreement is less good for the data pOl' which 
are distinctly smaller and start to increase strongly at somewhat smaller temperatures, 
for T > 200 K. In case of 02, our data agree for T > 150 K qualitatively well with the 
set [|0], correctly reproducing the gradual decrease with temperature, although our results 



are somewhat smaller. The data |19] exhibit here a different behavior, being of the same 



magnitude as the ones in Ref. |20|, but markedly fiat in the whole range of temperatures. 



We note here that in the work ||2T[, 0:2 < has been experimentally observed just below the 
melting point. On the other hand, the classical result for 02 calculated in Ref. exhibits 
instead an upward curvature. Concerning 0^3, our results agree quantitatively well with the 
data in the whole range of temperatures, although some scatter of the data precludes 
here a more detailed comparison. In the set [PD|, as behaves for T > 100 K qualitatively 
similar to our results, however, the plateau value is somewhat smaller and the strong increase 
in magnitude appears at a lower temperature, already between 200 - 250 K. The origin of 
some of the observed discrepancies may be either in the force field itself, or in the quantum 
effects, which may be relevant even at temperatures as high as 300 K 0. Also our fit of 



the data |2y] is not unique, and may itself be a source of additional errors in the coefficients 
a,. Last, but not least, it appears that there exists also a considerable scatter between the 
experimental data from different sources. We therefore believe that it would be interesting 
to re-examine the temperature dependence of the structural parameters (including possibly 
the setting angle of the chains) of well-crystalline samples of PE, in the whole range from 
very low temperatures up to the melting point, using up-to-date X-ray or neutron diffraction 
techniques. 
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In Fig. 12, the average fluctuations y of the setting angle of the chains are shown 

as a function of temperature. Apart from a trivial finite-size effect of average fluctuation 
increasing with decreasing chain length, we note that for the system with C12 chains the 
curve exhibits a marked enhancement of the fluctuations at T = 300 K. Such enhancement 
suggests that the system is approaching a transition into a phase similar to the "rotator" 
phases of n-parafflns, where the setting angle of the chains jumps among several minima. 
For the system with Cge chains, the same phenomenon occurs only for T > 400 K, which 
in a real PE crystal coincides with the melting point. In Fig. 13, we show a histogram of 
the setting angle ip for the system with Cge chains at the highest temperature T — 450 K. 
The distribution is still bimodal, with two peaks centered at id-^l), as it is in the ground 
state, corresponding to the "herringbone" arrangement of the chains. This proves that in our 
super-cell with periodic boundary conditions in all directions, the superheated orthorhombic 
structure is stable even at such high temperature, at least on the MC time scale of our run. 

We comment now briefly on some other internal structural parameters. In Figs. 14 and 
15, we show the average angles {6ccc) and {9hch) as a function of temperature. While the 
angle {Occc) develops between 250 and 300 K a very shallow minimum, the angle {Ohch) 
decreases roughly linearly throughout the whole range of temperatures, its overall variation 
being about a factor of 3 larger than that of the former. Both average bond lengths (rcc) 
and {tch) increase very slightly with temperature, the dependence being very close to linear, 
(rcc) varying from 1.5357 A at 10 K to 1.5398 A at 450 K, and {tch) varying from 1.0898 
A at 10 K to 1.0924 A at 450 K. It is interesting to show the average torsional fluctuations 
"cccc) ^ ^ function of temperature. Fig. 16, where we see an enhancement for T > 400 
K, corresponding to already discussed softening of the effective torsional potential (Fig. 11). 
The typical ffuctuation in this region is about 13°. From the histograms of the torsional 
angles we found that for temperatures up to 400 K, no gauche defects are created in the 
chains, while at 450 K only an extremely low population starts to arise. This is clearly a 
consequence of the periodic boundary conditions used in the chain direction together with 
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the still relatively short length of the chains used. 

We come now to some thermodynamic parameters, and start with the specific heat per 
unit cell. This is shown in Fig. 17 as a function of temperature. At low temperatures, where 
the classical crystal is always harmonic, ^ reaches the classical equipartition value of 18, 
corresponding to 3 x 12 = 36 degrees of freedom per unit cell. With increasing temperature, 
as the anharmonicities become important, its value gradually increases, until a finite size 
effect starts to appear for T > 200 K. While for the system with C12 chains Cp starts to 
grow faster for T > 200 K, for the system with Cge chains this occurs only at about 300 K. 
This behavior is likely to be related to the already discussed chain-length-dependent onset 
of increased fiuctuations of the setting angle, or rotations of the chains around their axes. 

In Figs. 18 - 22, we show the elastic constants Cn, C22, C12, C13, C23, determined by means 
of the Parrinello- Rahman fiuctuation formula (|l|), as a function of temperature. For the 
case of C33, we show in Fig. 23 both the values obtained according to the Parrinello- Rahman 
fiuctuation formula (|1|) and those found from the new formula (^. We see that the values 
obtained with the latter one have considerably smaller scatter. Actually, we have tried to 
use the formula (0) also for other elastic constants, however, only in the case of C33 a definite 
improvement of the convergence was observed. It would certainly be desirable to be able to 
improve on the accuracy of the elastic constants, in particular in case of 013,023; however, 
at present this seems to require a prohibitively large CPU time unless a considerably more 
efficient algorithm is available. At this point we note that some of the error bars shown in 
the figures for the elastic constants are probably an underestimate of the true ones, as they 
have been obtained from the variance of the results from just four batches. 

To discuss now the temperature dependence of the elastic constants, it's again convenient 
to do so separately for the lateral ones Cn, C22, C12, and for the ones related to the axial strain 
€3, namely C13, C23, C33. All three constants Cn, C22, C12 are seen to monotonously decrease with 
temperature, reaching at 400 K, just below the melting point, about 60 % or even less of 
their ground state values. This behavior originates mainly from the thermal expansion of 
the crystal in the lateral directions and is typical for van der Waals systems, like, e.g., solid 
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argon [|22|. An interesting finite-size effect is seen on the diagonal elastic constant C22, which 
is at T = 300 K distinctly smaller for the system with C12 chains than for other systems. 
This observation correlates with the effect seen on lattice parameter b, which was in turn 
found to be larger for the system with the shortest chains. 

Concerning the other group of elastic constants, we first note that the diagonal one C33 
is two orders of magnitude larger than all the other elastic constants. At low temperatures, 
C33 reaches a limiting value as large as 340 GPa and decreases monotonously with increasing 
temperature, dropping at 400 K to about 80 % of its ground state value. On the other hand, 
the two off-diagonal elastic constants, C13 and C23, are, interestingly, found to increase with 
temperature, in agreement with the results in Ref. 0. Taking into account the negative 
thermal expansion of the system in the axial direction, the behavior of this group of elastic 
constants does not appear to be just a trivial consequence of the thermal expansion, as it 
was in the case of cn, C22, cu. In order to understand such behavior, it is interesting to plot 
also the elastic constant C33 vs. the mean-square fluctuation of the torsional angle {<Pcccc)^ 
Fig. 24. We find again a roughly linear dependence (compare Fig. 10), which suggests that the 
softening of C33 is directly related to the activation of the torsional degrees of freedom. The 
mechanism underlying the behavior of these elastic constants then might be the following. 
At T = 0, where the chain backbones are perfectly fiat in their all-trans states, an axial 
deformation can be accommodated only by bending the bond angles Occc- Angle bending 
is, however, after bond stretching, the second most stiff degree of freedom in the system 
and determines the high ground state value of C33. As the torsional fluctuations become 
activated with increasing temperature, the chains start to "wiggle" and develop transverse 
fluctuations (see Fig.l (b)). In such configurations, it becomes possible to accommodate a 
part of the axial deformation in the torsional modes, which represent the most soft ones in 
the bonded interaction, and C33 is therefore renormalized to a smaller value. This must be, 
however, accompanied by an increased response of the transverse fluctuations of the chains, 
resulting in lateral strain, because of the "wiggling", since the total length of the chains is 
very hard to change. An increased lateral strain response to an axial strain means, however, 
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just an increase of the value of the elastic constants C13 and C23. It would be, of course, 
desirable to have a quantitative theory supporting this intuitive, but, as we believe, well 
plausible interpretation. 

The reason why the constants C13 and C23 appear to have the largest error bars among 
all elastic constants also seems to be connected with the fact that these two ones express 
just the coupling between the lateral and axial strains. The computational efficiency of the 
algorithm in determination of these two quantities is crucially dependent on the exchange of 
energy between the non-bonded and bonded interactions, which is probably still the most 
difficult point even with the present algorithm, because of the large separation of the relevant 
energy scales. 

Before closing this section, we would like to make yet few more remarks on the finite-size 
effects. Comparing the results at T = 300 K for both systems consisting of C24 chains, 
containing respectively 2 x 3 x 12 and 4 x 6 x 12 unit cells, we see that the values obtained 
with both system sizes are for all quantities practically equal. This suggests that the finite- 
size effects related directly to the volume of the box are relatively small, at least for the 
two system sizes considered (which are, however, still not very large). On the other hand, a 
definite chain-length-dependent finite-size effect has been found in case of several quantities 
for the chain lengths considered, its magnitude being also distinctly temperature dependent. 
For practical purposes, it results that the smallest system size used with C12 chains can be 
representative of a classical PE crystal only at rather low temperatures, perhaps below 
100 K, since at and above this temperature it exhibits pronounced finite-size effects. On the 
other hand, the systems with C24 and C48 chains appear to represent the classical PE crystal 
reasonably well for temperatures lower than 300 K, while at this and higher temperatures 
the use of a system with Cqq chains or even longer would be strongly recommended. 
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IV. CONCLUSIONS 



In this paper, we have demonstrated three main points. First, the MC algorithm using 
global moves on the chains in addition to the local moves on the atoms is a well applicable 
method for a classical simulation of crystalline PE. It allows an accurate determination of the 
structural properties and yields also fairly accurate results for the elastic constants. Second, 
the force field we have used [|TB[ is well able to reproduce the experimentally known structure 
in the whole range of temperatures, and where the classical description is appropriate, the 
results obtained agree well with the available experimental data. Third, we have studied 
the finite-size effects, mainly due to chain length, and determined for different temperatures 
a minimal chain length necessary for the system to be in the limit of long chains, and thus 
representative of PE. All these findings look promising for further studies, and here we 
suggest some possible directions. 

Basically, there are two routes to extend the present study, concentrating on the low- 
temperature and high-temperature region, respectively. The first one would aim on taking 
into account the quantum effects, e.g. by means of a Path Integral MC technique. Apart 
from improving the agreement with experiment, mainly at low temperatures, this technique 
is able to treat the quantum effects at a finite temperature in an essentially exact way and 
therefore should also allow to check the range of validity of various approximate treatments, 
like the quasi-harmonic or self-consistent quasi-harmonic approximations [0,^. This would 
help to understand better the true importance of quantum effects at different temperatures 
in this paradigmatic crystalline polymer system. 

The second route would aim on the high-temperature region of the phase diagram, where 
the orthorhombic crystal melts under normal pressure, but is known to undergo a phase 
transition into a hexagonal "condis" phase at elevated pressure 0. In the present simulation 
arrangement, melting is prevented by periodic boundary conditions in both lateral directions, 
but the same boundary conditions being applied in the chain direction inhibit also the 
creation of conformational defects. Nevertheless, there are indications arising from the 
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present simulation that a similar transition could indeed occur in the temperature region 
500 - 550 K, the main ones being the approach of the aspect ratio | towards the hexagonal 
value of \/3, and enhancement of the torsional angle and setting angle fluctuations for 
temperatures over 400 K. In order to study this high-temperature regime properly, several 
modifications of the present algorithm would be necessary. The main one would be a lifting 
of the periodic boundary conditions in the chain direction, thus introducing free chain ends. 
A use of a larger cutoff for the non-bonded interactions would be necessary in order to 
treat correctly the large amplitude displacements and rotations of the chains involved in the 
transition, and perhaps a non-spherical cutoff including certain number of atoms from one 
or two neighboring shells of chains around a given atom might be a preferred solution. It 
would also be necessary to modify the torsional potential in order to yield a correct value 
for the energy difference between the trans and gauche states. In order to be able to sample 
configurations with a considerable population of conformational defects, it might also be 
useful to introduce some other kind of MC moves acting on the torsional degrees of freedom. 
Finally, a full MC version of the Parrinello- Rahman variable-cell-shape technique should 
be introduced in order to allow also for shear fiuctuations. These might in principle be 
substantially involved in the transition itself, although the hexagonal phase can be reached 
from the orthorhombic one without creating a static shear strain. 

Before closing, we also emphasize that techniques similar to those applied here should 
be useful for a wide variety of other macromolecular crystals. 
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FIGURES 



Fig.l. A snapshot of a projection of a typical configuration of the system at T = 400 
K on the xy plane (a), and on the xz plane (b). The size of the reference super-cell is 

= 14. SA, Ly = 15. OA, Lz = 121. 44A. Filled points represent the carbon atoms, empty 
ones the hydrogen atoms. Note that the "herringbone" arrangement of the chains is still 
clearly visible in (a). 



Fig. 2. Temperature dependence of the lattice parameter a, shown for different system 
sizes, together with the experimental data |[T9|j20[] . In this and most of the following figures, 
statistical error bars are shown, lines are for visual help only. 



Fig. 3. Temperature dependence of the lattice parameter 6, shown for different system 



sizes, together with the experimental data ||I9|,E0|. Note the finite-size effect 



Fig. 4. Temperature dependence of the lattice parameter c, shown for different system 
sizes, together with the experimental data 1|T^^ . 



Fig. 5. Temperature dependence of the aspect ratio |. Note the approach towards the 
hexagonal value of = 1.732 with increasing temperature. 

Fig. 6. Temperature dependence of the average setting angle (|-?/'|) of the chains, shown 
for different system sizes. Note the finite-size effect. 



Fig. 7. Thermal expansion coefficient ai as a function of temperature, shown for different 
system sizes, together with the experimental data [|1^,^ . 
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Fig. 8. Thermal expansion coefficient 02 as a function of temperature, shown for different 
system sizes, together with the experimental data [0,^]. Note the finite-size effect. 



Fig. 9. Thermal expansion coefficient as as a function of temperature, shown for different 



Fig. 10. Lattice parameter c vs. {<f>c'ccc)y shown for different system sizes. Note the 
linear dependence over the whole temperature range, which can be fitted by a line c = 



Fig. 11. Temperature dependence of the effective torsional force constant Ctors = 



shown for different system sizes. Note the considerable softening with increasing tempera- 
ture. 



Fig. 12. Temperature dependence of the average fluctuation \J {{S\ip\y) of the setting 
angle of the chains, shown for different system sizes. 

Fig. 13. Histogram of the setting angle ip of the chains for the system with Cge chains at 
the highest temperature T = 450 K. Note that the distribution is still bimodal, corresponding 
to the "herringbone" arrangement of the chains. 

Fig. 14. Temperature dependence of the average bond angle {Occc)i shown for different 
system sizes. Note the minimum at 300 K. 

Fig. 15. Temperature dependence of the average bond angle {9hch), shown for different 
system sizes. 



system sizes, together with the experimental data [0,^]. 



2.530- 9.50 X 10~H(i)r 




cccc 
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Fig. 16. Temperature dependence of the average fluctuation \l {(fcccc) torsional 
angle ^cccc from the trans minimum, shown for different system sizes. 

Fig. 17. Temperature dependence of the speciflc heat per unit cell at constant pressure 
Cp, shown for different system sizes. 

Fig. 18. Elastic constant Cn as a function of temperature, shown for different system 
sizes. 

Fig. 19. Elastic constant C22 as a function of temperature, shown for different system 
sizes. Note the finite-size effect. 

Fig. 20. Elastic constant C12 as a function of temperature, shown for different system 
sizes. 

Fig. 21. Elastic constant C13 as a function of temperature, shown for different system 
sizes. Note the increase with temperature. 

Fig. 22. Elastic constant C23 as a function of temperature, shown for different system 
sizes. Note the increase with temperature. 

Fig. 23. Elastic constant C33 as a function of temperature, shown for different system 
sizes. Empty and filled symbols represent data determined from the formulas (|I|) and @, 
respectively. Note the smaller scatter of the latter results. 

Fig. 24. Elastic constant C33 vs. {(j^cccc)^ shown for different system sizes. Note the 
roughly linear dependence. 
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